/********************************************************************/
/********************************************************************/
/** COMPUTE COMMAND LINE SCIENTIFIC CALCULATOR                     **/
/********************************************************************/
/********************************************************************/
/**                                                                **/
/**       IBM Corporation                                          **/
/**       HARRISMH at RCHVMP2                                      **/
/**       EL8/663-1 B627                                           **/
/**       Rochester, MN  55901                                     **/
/**                                                                **/
/**       AUTHOR:  Mark H. Harris                                  **/
/**         DATE:  93/10/22                                        **/
/**       UPDATE:  98/03/07                                        **/
/**          VER:  2.0a                                            **/
/**                                                                **/
/********************************************************************/
/********************************************************************/
/**                          syntax                                **/
/**                                                                **/
/**  COMPUTE  expression<;expression><;expression><@digits>        **/
/**                                                                **/
/**                                                                **/
/**  The expression(s) will be computed and displayed in terminal  **/
/**  line mode to arbitrary <@digits> of accuracy.                 **/
/**                                                                **/
/**  The expression(s) may be any REXX math clause and may include **/
/**  a variable assignment, for use in a subsequent expression:    **/
/**  ie.,  COMPUTE A=3;B=7;A/B;A*B@30                              **/
/**        (will compute both the quotient and product of A & B)   **/
/**                                                                **/
/********************************************************************/
/********************************************************************/
/**                          NOTES                                 **/
/**                                                                **/
/**  COMPUTE  expression<;expression><;expression><@digits>        **/
/**                                                                **/
/**  This program exploits Rexx  INTERPRET and NUMERIC DIGITS      **/
/**                                                                **/
/**  This exec is portable to OS/2, NT and w95 Rexx/objectRexx     **/
/**  without changes.                                              **/
/**                                                                **/
/********************************************************************/
/********************************************************************/
/**  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND        **/
/**  CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,   **/
/**  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF      **/
/**  MECHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE       **/
/**  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR         **/
/**  CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,  **/
/**  SPECIAL, EXEMPLARY, OR CONSEQUENCIAL DAMAGES (INCLUDING, BUT  **/
/**  NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;  **/
/**  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERUPTION)       **/
/**  HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN     **/
/**  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR  **/
/**  OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE **/
/**  EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.            **/
/********************************************************************/
/********************************************************************/
/**                          updates                               **/
/** 93/10/21 mhh New Program                                       **/
/** 98/03/07 mhh New Version 2.0a                                  **/
/**              This version includes trigonometric, logarithmic, **/
/**              exponential and combinatorial functions.          **/
/**                                                                **/
/**              See the REXROOTS ABOUT file for a discussion      **/
/**              of the included functions and some examples,      **/
/**              using the COMPUTE program.                        **/
/**                                                                **/
/********************************************************************/
/********************************************************************/
/**                          dependencies                          **/
/**                                                                **/
/** Numeric Fuzz:                                                  **/
/**              Most of these routines bump-up the numeric digits **/
/**              by a fuzzFactor.  Rexx fuzz is used to terminate  **/
/**              the do-forever iteration at the correct number of **/
/**              required places.                                  **/
/**              The routine then sets the numeric digits back to  **/
/**              the application setting and returns result-0.     **/
/**              The (result-0) forces truncation/rounding to give **/
/**              the expected result.                              **/
/**              Some implementations of REXX do not support the   **/
/**              fuzz option;  the iterative do-forever methods    **/
/**              would need modification for correct termination.  **/
/**                                                                **/
/********************************************************************/
/********************************************************************/
/**                          mainline                              **/
/********************************************************************/
arg expression '@' max_digits .
   select
      when expression='' | expression='?' then say syntax_diagram()
      when datatype(max_digits,'W') & max_digits>0 then,
                                          numeric digits max_digits
      otherwise
         numeric digits 18
      end
      call evaluate expression
   exit 0
/********************************************************************/
/**                          procedures                            **/
/********************************************************************/
SYNTAX_DIAGRAM: procedure
   return'Syntax:   COMPUTE expression <;expression><@precision>'
 
EVALUATE: procedure expose max_digits
arg expression
signal on error; signal on syntax
   /* 502 digits pi/2 */
   hpi="1.57079632679489661923132169163975144209858469968755291048747229615390820314310449931401741267105853"
   hpi=hpi"3991074043256641153323546922304775291115862679704064240558725142051350969260552779822311474477465190"
   hpi=hpi"9822144054878329667230642378241168933915826356009545728242834617301743052271633241066968036301245706"
   hpi=hpi"3686229350330315779408744076046048141462704585768218394629518000566526527441023326069207347597075580"
   hpi=hpi"471652863518287979597654609305869096630589655255927403723118998137478367594287636244561396909150597456"
   /* 502 digits pi */
   pi="3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679"
   pi=pi"8214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196"
   pi=pi"4428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273"
   pi=pi"7245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094"
   pi=pi"3305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912"
 
   resetTime=time('R')
   do while expression <> ''
      parse value expression with exp1';'expression
      select
         when pos('=',exp1)<>0 then do
            interpret 'do; 'exp1'; end;'
            end
         otherwise
            interpret 'do; answer='exp1'; end;'
            select
               when answer='' then say syntax_diagram()
               when datatype(answer,'N') then say answer
               otherwise signal error
               end
         end
      end
   say "elapsed time:" time('E')
   return
/********************************************************************/
/**                          routines                              **/
/********************************************************************/
ERROR:; SYNTAX:
   say 'Does Not Compute'
   exit 28
/********************************************************************/
/**                          functions                             **/
/********************************************************************/
 
fx: procedure /* template
                                                 */
arg x
   maxDigits=digits()
   numeric digits maxDigits*2
   numeric fuzz maxDigits
      /* f(x) goes here */
   numeric fuzz 0
   numeric digits maxDigits
   return x-0
 
/********************************************************************/
/** Exponential - Logarithmic                                      **/
/********************************************************************/
exp: procedure /* (a>=0) (n>=0)
                                         */
arg a,n
   select
      when datatype(n, "W")=1 then x=a**n
      otherwise x=e(ln(a)*n)
      end
   return x
 
 
ln: procedure /* series methods ( x > 0 ) */
arg x
   if x<=0 then return "ERROR"
   maxDigits=digits()
   numeric digits maxDigits+7
   numeric fuzz 7
   select
      when x>80 then do
         i=0
         do while x>=5
            i=i+1
            x=x/10
            end
         x=ln(x)+i*ln(10)
         end
      when x<.5 then do
         i=0
         do while x<.5
            i=i+1
            x=x*10
            end
         x=ln(x)-i*ln(10)
         end
      when x>1.6 then do
         u=(x-1)/(x+1); x=u; c=1; s=u**2
         do forever
            c=c+1
            u=s*u
            x1=x+u/(2*c-1)
            if x1=x then leave
            x=x1
            end
         x=x*2
         end
      otherwise
         u=(x-1); x=u; c=1; s=u; sign=1
         do forever
            c=c+1
            u=u*s
            sign=sign*(-1)
            x1=x+sign*u/c
            if x1=x then leave
            x=x1
            end
      end
   numeric fuzz 0
   numeric digits maxDigits
   return x-0
 
e: procedure /* (all real x)
   e^x = 1+x+x^2/2!+x^3/3!-x^4/4!... */
arg x
   maxDigits=digits()
   select
      when x<(-5) then do
         fuzzFactor=abs(x)%5
         numeric digits maxDigits*fuzzFactor
         numeric fuzz maxDigits*(fuzzFactor-1)
         end
      otherwise
         numeric digits maxDigits+7
         numeric fuzz 7
      end
   n=x; s=x; x=x+1; d=1; c=1
   do forever
      c=c+1
      n=n*s
      d=d*c
      x1=x+n/d
      if x1=x then leave
      x=x1
      end
   numeric fuzz 0
   numeric digits maxDigits
   return x-0
 
/********************************************************************/
/** Trigonometric                                                  **/
/********************************************************************/
tan: procedure /* (all real x)
                                      */
arg x,hyper
   maxDigits=digits()
   numeric digits maxDigits*2
   numeric fuzz maxDigits
   if hyper="" then x=sin(x)/cos(x)
   else x=sinh(x)/cosh(x)
   numeric fuzz 0
   numeric digits maxDigits
   return x-0
 
tanh: procedure /* (all real x)
                                    */
arg x
   return tan(x,"h")
 
cos: procedure /* (all real x)
     cos(x)=1-x^2/2!+x^4/4!-x^6/6!...
    cosh(x)=1+x^2/2!+x^4/4!+x^6/6!... */
 
arg x,hyper
   maxDigits=digits()
   numeric digits maxDigits*2
   numeric fuzz maxDigits
   n=1; s=x**2; x=1; sign=1; d=1; c=1
   select
      when hyper="" then do forever      /* cosine function */
         c=c+1
         n=n*s
         sign=(-1)*sign
         do i=2 to 3 by 1
            d=d*(2*c-i)
            end
         x1=x+sign*n/d
         if x1=x then leave
         x=x1
         end
      otherwise
         do forever       /* hyperbolic cosine function */
            c=c+1
            n=n*s
            do i=2 to 3 by 1
               d=d*(2*c-i)
               end
            x1=x+n/d
            if x1=x then leave
            x=x1
            end
      end
   numeric fuzz 0
   numeric digits maxDigits
   return x-0
 
cosh: procedure /* (all real x)
                                           */
arg x
   return cos(x,"h")
 
sin: procedure /* (all real x)
     sin(x)=x/1!-x^3/3!+x^5/5!-x^7/7!...
    sinh(x)=x/1!+x^3/3!+x^5/5!+x^7/7!...   */
 
arg x,hyper
   maxDigits=digits()
   numeric digits maxDigits*2
   numeric fuzz maxDigits
   n=x; s=x**2; sign=1; d=1; c=1
   select
      when hyper="" then do forever       /* sine function */
         c=c+1
         n=n*s
         sign=(-1)*sign
         do i=1 to 2 by 1
            d=d*(2*c-i)
            end
         x1=x+sign*n/d
         if X1=x then leave
         x=x1
         end
      otherwise
         do forever      /* hyperbolic sine Function */
            c=c+1
            n=n*s
            do i=1 to 2 by 1
               d=d*(2*c-i)
               end
            x1=x+n/d
            if X1=x then leave
            x=x1
            end
      end /*select*/
   numeric fuzz 0
   numeric digits maxDigits
   return x-0
 
sinh: procedure /* (all real x)
                                            */
arg x
   return sin(x,"h")
 
sin_1: procedure expose hpi /* ( x^2<1  -pi/2 < sin_1(x) < pi/2 )
       sin_1(x)=x+(x^3/2*3)+(x^5*1*3/2*4*5)+(x^7*1*3*5/2*4*6*7)... */
arg x
   if x**2 >1 then return "ERROR"
   maxDigits=digits()
   numeric digits maxDigits*2
   numeric fuzz maxDigits
   n=x; d=1; s=x**2; c=1; p=0
   if abs(x)>(.7854) then do
         p=1; if x<0 then p=-1
         s=1-s; x=sqrt(s); n=x
         end
   do forever
      c=c+1; cc=2*c
      n=n*(cc-3)*s
      d=d*(cc-2)
      x1=x+n/(d*(cc-1))
      if x1=x then leave
      x=x1
      end
   numeric fuzz 0
   numeric digits maxDigits
   select
      when p>0 then return hpi-x
      when p<0 then return x-hpi
      otherwise nop
      end
   return x-0
 
cos_1: procedure expose hpi /*  ( x^2<1  0 < cos_1(x) < pi )
                                      */
arg x
   return hpi-sin_1(x)
 
tan_1: procedure expose hpi /*  ( all real x )
                                      */
arg x
   return hpi-cos_1(x/sqrt(x**2+1))
 
/********************************************************************/
/** Combinatorial                                                  **/
/********************************************************************/
fact: procedure /* n!  n=(0,1,2,3...) */
arg n
   numeric digits 500
   m=1
   do i=2 to n by 1
      m=i*m
      end
   return m
 
part: procedure /* Ordered Partition */
arg m parts
   numeric digits 500
   mFact=fact(m)
   do while parts<>''
      parse var parts part parts
      if datatype(part,'W') & m-part>=0 then do
         m=m-part
         mFact=mFact/fact(part)
         end
      end
   return mFact/fact(m)
 
/********************************************************************/
/**  PI routines                                                   **/
/********************************************************************/
PI: procedure /* AGM Gauss-Legendre Method
    Crandall:  Projects In Scientific Computation */
arg n
   maxDigits=digits()
   numeric digits maxDigits+7
   numeric fuzz 7
   x=sqrt(2); p=2+x; y=sqrt(x)
   do forever
      s=sqrt(x)
      x=(s+1/s)/2
      p1=p*((x+1)/(y+1))
      if p1=p then leave
      p=p1
      s=sqrt(x)
      y=((y*s)+(1/s))/(y+1)
      end
   numeric fuzz 0
   numeric digits maxDigits
   return p-0
 
PIa: procedure /* AGM Gauss-Legendre Method */
arg n
   maxDigits=digits()
   numeric digits maxDigits+7
   numeric fuzz 7
   a=1; b=1/sqrt(2); t=1/4; x=1
   do i=1 to 100
      y=a
      a=(a+b)/2
      b=sqrt(b*y)
      if a=b then leave
      t=t-x*(y-a)**2
      x=2*x
      end
   PI=((a+b)**2)/(4*t)
   numeric fuzz 0
   numeric digits maxDigits
   return PI-0
 
/********************************************************************/
/** Roots                                                          **/
/********************************************************************/
sqrt: procedure /* (all real x) */
arg x
   return nrt(x,2)
 
cbrt: procedure /* (all real x) */
arg x
   return nrt(x,3)
 
nrt: procedure /* ( a>=0 ) ( n>=0 )
                                     */
arg a,n
   maxDigits=digits()
   numeric digits maxDigits+7
   numeric fuzz 7
   if a=1 | a=0 then return a
   if n=0 then return 1
   if a<0 | n<0 then return "ERROR"
   select
      when datatype(n, "W" )=1 then do
         m=n-1
         x=bisect(a,n)
         do forever  /* Newton Method  x:= x - ( f(x)/f'(x) ) */
            x1=(m*x**n+a)/(n*x**m)
            if x1=x then leave
            x=x1
            end
         end
      otherwise
         x=e(ln(a)/n)
      end
   numeric fuzz 0
   numeric digits maxDigits
   return x-0
 
bisect: procedure /* Bisect Method used to seed nrt() */
arg n,root /* (n>=0) (root: 2,3,4,5,6...) */
   numeric digits 12
   numeric fuzz 3
   a=0; b=n
   if n<1 then b=1
   else a=1
   do forever
      m=(a+b)/2
      Fm=m**root-n
      select
         when Fm>0 then b=m
         when Fm<0 then a=m
         otherwise
            a=b
         end
      if a=b then leave
      end
   return b
